
R version 4.1.2 (2021-11-01) -- "Bird Hippie"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> ## Section 3.2.4.B of the pre-analysis plan: "Entropy of Posterior
> ## Probabilities of Ranks (conditional on coding task and data-collection
> ## procedure)"
> ##
> ## Kevin Quinn
> ## University of Michigan
> ##
> ## 7/15/2019
> ##
> 
> #library(TopKLists)
> library(coda)
> 
> set.seed(2102983)
> 
> 
> ########################################################################
> ## START a1-b1 comparisons
> ########################################################################
> 
> ## load MCMC output
> load("MM3.2.a1.M1.out.Rda")
> a1.M1.out <- a1.M1.out[sample(1:nrow(a1.M1.out), nrow(a1.M1.out),
+                               replace=FALSE), ]
> load("MM3.2.a1.M2.out.Rda")
> a1.M2.out <- a1.M2.out[sample(1:nrow(a1.M2.out), nrow(a1.M2.out),
+                               replace=FALSE), ]
> load("MM3.2.a1.pairwise.out.Rda")
> a1.pairwise.out <- a1.pairwise.out[sample(1:nrow(a1.pairwise.out),
+                                           nrow(a1.pairwise.out),
+                               replace=FALSE), ]
> 
> load("MM3.2.b1.M1.out.Rda")
> b1.M1.out <- b1.M1.out[sample(1:nrow(b1.M1.out), nrow(b1.M1.out),
+                               replace=FALSE), ]
> 
> load("MM3.2.b1.M2.out.Rda")
> b1.M2.out <- b1.M2.out[sample(1:nrow(b1.M2.out), nrow(b1.M2.out),
+                               replace=FALSE), ]
> load("MM3.2.b1.pairwise.out.Rda")
> b1.pairwise.out <- b1.pairwise.out[sample(1:nrow(b1.pairwise.out),
+                                           nrow(b1.pairwise.out),
+                               replace=FALSE), ]
> 
> ## make MCMC output comparable
> a1.M1.out <- a1.M1.out[,1:48]
> a1.M2.out <- a1.M2.out[,1:48]
> b1.M1.out <- b1.M1.out[,1:48]
> b1.M2.out <- b1.M2.out[,1:48]
> 
> colnames(a1.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a1.M1.out))
> colnames(a1.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a1.M2.out))
> colnames(b1.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b1.M1.out))
> colnames(b1.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b1.M2.out))
> colnames(a1.pairwise.out) <- gsub("theta\\.", "", colnames(a1.pairwise.out))
> colnames(b1.pairwise.out) <- gsub("theta\\.", "", colnames(b1.pairwise.out))
> 
> for (j in 2:48){
+     a1.M1.out[,j] <- a1.M1.out[,1] + a1.M1.out[,j]
+     a1.M2.out[,j] <- a1.M2.out[,1] + a1.M2.out[,j]
+     b1.M1.out[,j] <- b1.M1.out[,1] + b1.M1.out[,j]
+     b1.M2.out[,j] <- b1.M2.out[,1] + b1.M2.out[,j]
+ }
> colnames(a1.M1.out)[1] <- colnames(a1.pairwise.out)[1]
> colnames(a1.M2.out)[1] <- colnames(a1.pairwise.out)[1]
> colnames(b1.M1.out)[1] <- colnames(a1.pairwise.out)[1]
> colnames(b1.M2.out)[1] <- colnames(a1.pairwise.out)[1]
> 
> 
> M <- nrow(a1.M1.out)
> J <- ncol(a1.M1.out)
> 
> rank.a1.M1 <- matrix(NA, M, J)
> rank.a1.M2 <- matrix(NA, M, J)
> rank.b1.M1 <- matrix(NA, M, J)
> rank.b1.M2 <- matrix(NA, M, J)
> rank.a1.pair <- matrix(NA, M, J)
> rank.b1.pair <- matrix(NA, M, J)
> 
> for (iter in 1:M){
+     rank.a1.M1[iter,] <- rank(a1.M1.out[iter,])
+     rank.a1.M2[iter,] <- rank(a1.M2.out[iter,])
+     rank.b1.M1[iter,] <- rank(b1.M1.out[iter,])
+     rank.b1.M2[iter,] <- rank(b1.M2.out[iter,])
+     rank.a1.pair[iter,] <- rank(a1.pairwise.out[iter,])
+     rank.b1.pair[iter,] <- rank(b1.pairwise.out[iter,])
+ }
> 
> rank.probs.a1.M1 <- matrix(NA, J, J)
> rank.probs.a1.M2 <- matrix(NA, J, J)
> rank.probs.b1.M1 <- matrix(NA, J, J)
> rank.probs.b1.M2 <- matrix(NA, J, J)
> rank.probs.a1.pair <- matrix(NA, J, J)
> rank.probs.b1.pair <- matrix(NA, J, J)
> 
> for (j1 in 1:J){
+     for (j2 in 1:J){
+         rank.probs.a1.M1[j1,j2] <- mean(rank.a1.M1[, j1] == j2)
+         rank.probs.a1.M2[j1,j2] <- mean(rank.a1.M2[, j1] == j2)
+         rank.probs.b1.M1[j1,j2] <- mean(rank.b1.M1[, j1] == j2)
+         rank.probs.b1.M2[j1,j2] <- mean(rank.b1.M2[, j1] == j2)
+         rank.probs.a1.pair[j1,j2] <- mean(rank.a1.pair[, j1] == j2)
+         rank.probs.b1.pair[j1,j2] <- mean(rank.b1.pair[, j1] == j2)
+     }
+ }
> 
> h.a1.M1 <- rep(NA, J)
> h.a1.M2 <- rep(NA, J)
> h.b1.M1 <- rep(NA, J)
> h.b1.M2 <- rep(NA, J)
> h.a1.pair <- rep(NA, J)
> h.b1.pair <- rep(NA, J)
> for (j in 1:J){
+     pvec <- rank.probs.a1.M1[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a1.M1[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.a1.M2[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a1.M2[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b1.M1[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b1.M1[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b1.M2[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b1.M2[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.a1.pair[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a1.pair[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b1.pair[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b1.pair[j] <- -1 * sum(holder)
+     
+ }
> 
> rank.a1.pair.hat <- apply(rank.a1.pair, 2, mean)
> ordvec <- order(rank.a1.pair.hat)
> 
> pdf(file="MM3-2-4-B-a1.pdf", height=6, width=7)
> plot(h.a1.M1[ordvec], xlab="Photo Index",
+      ylab="Entropy",
+      main="Coder Race Correlates with Photos Observed (Column a)", pch=0,
+      col="orangered1", ylim=c(0,5))
> abline(h=0)
> points(h.a1.M1[ordvec], pch=0, col="orangered1")
> points(h.a1.M2[ordvec], pch=1, col="dodgerblue")
> points(h.a1.pair[ordvec], pch=3, col="purple4")
> legend(x=0, y=5, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
+        col=c("orangered1", "dodgerblue", "purple4"))
> dev.off()
null device 
          1 
> 
> 
> 
> rank.b1.pair.hat <- apply(rank.b1.pair, 2, mean)
> ordvec <- order(rank.b1.pair.hat)
> 
> pdf(file="MM3-2-4-B-b1.pdf", height=6, width=7)
> plot(h.b1.M1[ordvec], xlab="Photo Index",
+      ylab="Entropy",
+      main="Coder Race Correlates with Photos Observed (Column b)", pch=0,
+      col="orangered1", ylim=c(0,5))
> abline(h=0)
> points(h.b1.M1[ordvec], pch=0, col="orangered1")
> points(h.b1.M2[ordvec], pch=1, col="dodgerblue")
> points(h.b1.pair[ordvec], pch=3, col="purple4")
> legend(x=0, y=5, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
+        col=c("orangered1", "dodgerblue", "purple4"))
> dev.off()
null device 
          1 
> 
> 
> 
> 
> 
> ########################################################################
> ## END a1-b1 comparisons
> ########################################################################
> 
> 
> 
> 
> 
> 
> ########################################################################
> ## START a2-b2 comparisons
> ########################################################################
> 
> ## load MCMC output
> load("MM3.2.a2.M1.out.Rda")
> a2.M1.out <- a2.M1.out[sample(1:nrow(a2.M1.out), nrow(a2.M1.out),
+                               replace=FALSE), ]
> load("MM3.2.a2.M2.out.Rda")
> a2.M2.out <- a2.M2.out[sample(1:nrow(a2.M2.out), nrow(a2.M2.out),
+                               replace=FALSE), ]
> load("MM3.2.a2.pairwise.out.Rda")
> a2.pairwise.out <- a2.pairwise.out[sample(1:nrow(a2.pairwise.out),
+                                           nrow(a2.pairwise.out),
+                               replace=FALSE), ]
> 
> load("MM3.2.b2.M1.out.Rda")
> b2.M1.out <- b2.M1.out[sample(1:nrow(b2.M1.out), nrow(b2.M1.out),
+                               replace=FALSE), ]
> 
> load("MM3.2.b2.M2.out.Rda")
> b2.M2.out <- b2.M2.out[sample(1:nrow(b2.M2.out), nrow(b2.M2.out),
+                               replace=FALSE), ]
> load("MM3.2.b2.pairwise.out.Rda")
> b2.pairwise.out <- b2.pairwise.out[sample(1:nrow(b2.pairwise.out),
+                                           nrow(b2.pairwise.out),
+                               replace=FALSE), ]
> 
> ## make MCMC output comparable
> a2.M1.out <- a2.M1.out[,1:48]
> a2.M2.out <- a2.M2.out[,1:48]
> b2.M1.out <- b2.M1.out[,1:48]
> b2.M2.out <- b2.M2.out[,1:48]
> 
> colnames(a2.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a2.M1.out))
> colnames(a2.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a2.M2.out))
> colnames(b2.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b2.M1.out))
> colnames(b2.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b2.M2.out))
> colnames(a2.pairwise.out) <- gsub("theta\\.", "", colnames(a2.pairwise.out))
> colnames(b2.pairwise.out) <- gsub("theta\\.", "", colnames(b2.pairwise.out))
> 
> for (j in 2:48){
+     a2.M1.out[,j] <- a2.M1.out[,1] + a2.M1.out[,j]
+     a2.M2.out[,j] <- a2.M2.out[,1] + a2.M2.out[,j]
+     b2.M1.out[,j] <- b2.M1.out[,1] + b2.M1.out[,j]
+     b2.M2.out[,j] <- b2.M2.out[,1] + b2.M2.out[,j]
+ }
> colnames(a2.M1.out)[1] <- colnames(a2.pairwise.out)[1]
> colnames(a2.M2.out)[1] <- colnames(a2.pairwise.out)[1]
> colnames(b2.M1.out)[1] <- colnames(a2.pairwise.out)[1]
> colnames(b2.M2.out)[1] <- colnames(a2.pairwise.out)[1]
> 
> 
> M <- nrow(a2.M1.out)
> J <- ncol(a1.M1.out)
> 
> rank.a2.M1 <- matrix(NA, M, J)
> rank.a2.M2 <- matrix(NA, M, J)
> rank.b2.M1 <- matrix(NA, M, J)
> rank.b2.M2 <- matrix(NA, M, J)
> rank.a2.pair <- matrix(NA, M, J)
> rank.b2.pair <- matrix(NA, M, J)
> 
> for (iter in 1:M){
+     rank.a2.M1[iter,] <- rank(a2.M1.out[iter,])
+     rank.a2.M2[iter,] <- rank(a2.M2.out[iter,])
+     rank.b2.M1[iter,] <- rank(b2.M1.out[iter,])
+     rank.b2.M2[iter,] <- rank(b2.M2.out[iter,])
+     rank.a2.pair[iter,] <- rank(a2.pairwise.out[iter,])
+     rank.b2.pair[iter,] <- rank(b2.pairwise.out[iter,])
+ }
> 
> 
> rank.probs.a2.M1 <- matrix(NA, J, J)
> rank.probs.a2.M2 <- matrix(NA, J, J)
> rank.probs.b2.M1 <- matrix(NA, J, J)
> rank.probs.b2.M2 <- matrix(NA, J, J)
> rank.probs.a2.pair <- matrix(NA, J, J)
> rank.probs.b2.pair <- matrix(NA, J, J)
> 
> for (j1 in 1:J){
+     for (j2 in 1:J){
+         rank.probs.a2.M1[j1,j2] <- mean(rank.a2.M1[, j1] == j2)
+         rank.probs.a2.M2[j1,j2] <- mean(rank.a2.M2[, j1] == j2)
+         rank.probs.b2.M1[j1,j2] <- mean(rank.b2.M1[, j1] == j2)
+         rank.probs.b2.M2[j1,j2] <- mean(rank.b2.M2[, j1] == j2)
+         rank.probs.a2.pair[j1,j2] <- mean(rank.a2.pair[, j1] == j2)
+         rank.probs.b2.pair[j1,j2] <- mean(rank.b2.pair[, j1] == j2)
+     }
+ }
> 
> h.a2.M1 <- rep(NA, J)
> h.a2.M2 <- rep(NA, J)
> h.b2.M1 <- rep(NA, J)
> h.b2.M2 <- rep(NA, J)
> h.a2.pair <- rep(NA, J)
> h.b2.pair <- rep(NA, J)
> for (j in 1:J){
+     pvec <- rank.probs.a2.M1[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a2.M1[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.a2.M2[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a2.M2[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b2.M1[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b2.M1[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b2.M2[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b2.M2[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.a2.pair[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a2.pair[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b2.pair[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b2.pair[j] <- -1 * sum(holder)
+     
+ }
> 
> rank.a2.pair.hat <- apply(rank.a2.pair, 2, mean)
> ordvec <- order(rank.a2.pair.hat)
> 
> pdf(file="MM3-2-4-B-a2.pdf", height=6, width=7)
> plot(h.a2.M1[ordvec], xlab="Photo Index",
+      ylab="Entropy",
+      main="Distractor-Photo Race Correlates with Photos Observed (Column a)", pch=0,
+      col="orangered1", ylim=c(0,5))
> abline(h=0)
> points(h.a2.M1[ordvec], pch=0, col="orangered1")
> points(h.a2.M2[ordvec], pch=1, col="dodgerblue")
> points(h.a2.pair[ordvec], pch=3, col="purple4")
> legend(x=0, y=5, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
+        col=c("orangered1", "dodgerblue", "purple4"))
> dev.off()
null device 
          1 
> 
> 
> 
> rank.b2.pair.hat <- apply(rank.b2.pair, 2, mean)
> ordvec <- order(rank.b2.pair.hat)
> 
> pdf(file="MM3-2-4-B-b2.pdf", height=6, width=7)
> plot(h.b2.M1[ordvec], xlab="Photo Index",
+      ylab="Entropy",
+      main="Distractor-Photo Race Correlates with Photos Observed (Column b)", pch=0,
+      col="orangered1", ylim=c(0,5))
> abline(h=0)
> points(h.b2.M1[ordvec], pch=0, col="orangered1")
> points(h.b2.M2[ordvec], pch=1, col="dodgerblue")
> points(h.b2.pair[ordvec], pch=3, col="purple4")
> legend(x=0, y=5, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
+        col=c("orangered1", "dodgerblue", "purple4"))
> dev.off()
null device 
          1 
> 
> 
> 
> ########################################################################
> ## END a2-b2 comparisons
> ########################################################################
> 
> 
> 
> 
> 
> 
> ########################################################################
> ## START a3-b3 comparisons
> ########################################################################
> 
> ## load MCMC output
> load("MM3.2.a3.M1.out.Rda")
> a3.M1.out <- a3.M1.out[sample(1:nrow(a3.M1.out), nrow(a3.M1.out),
+                               replace=FALSE), ]
> load("MM3.2.a3.M2.out.Rda")
> a3.M2.out <- a3.M2.out[sample(1:nrow(a3.M2.out), nrow(a3.M2.out),
+                               replace=FALSE), ]
> load("MM3.2.a3.pairwise.out.Rda")
> a3.pairwise.out <- a3.pairwise.out[sample(1:nrow(a3.pairwise.out),
+                                           nrow(a3.pairwise.out),
+                               replace=FALSE), ]
> 
> load("MM3.2.b3.M1.out.Rda")
> b3.M1.out <- b3.M1.out[sample(1:nrow(b3.M1.out), nrow(b3.M1.out),
+                               replace=FALSE), ]
> 
> load("MM3.2.b3.M2.out.Rda")
> b3.M2.out <- b3.M2.out[sample(1:nrow(b3.M2.out), nrow(b3.M2.out),
+                               replace=FALSE), ]
> load("MM3.2.b3.pairwise.out.Rda")
> b3.pairwise.out <- b3.pairwise.out[sample(1:nrow(b3.pairwise.out),
+                                           nrow(b3.pairwise.out),
+                               replace=FALSE), ]
> 
> ## make MCMC output comparable
> a3.M1.out <- a3.M1.out[,1:48]
> a3.M2.out <- a3.M2.out[,1:48]
> b3.M1.out <- b3.M1.out[,1:48]
> b3.M2.out <- b3.M2.out[,1:48]
> 
> colnames(a3.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a3.M1.out))
> colnames(a3.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a3.M2.out))
> colnames(b3.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b3.M1.out))
> colnames(b3.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b3.M2.out))
> colnames(a3.pairwise.out) <- gsub("theta\\.", "", colnames(a3.pairwise.out))
> colnames(b3.pairwise.out) <- gsub("theta\\.", "", colnames(b3.pairwise.out))
> 
> for (j in 2:48){
+     a3.M1.out[,j] <- a3.M1.out[,1] + a3.M1.out[,j]
+     a3.M2.out[,j] <- a3.M2.out[,1] + a3.M2.out[,j]
+     b3.M1.out[,j] <- b3.M1.out[,1] + b3.M1.out[,j]
+     b3.M2.out[,j] <- b3.M2.out[,1] + b3.M2.out[,j]
+ }
> colnames(a3.M1.out)[1] <- colnames(a3.pairwise.out)[1]
> colnames(a3.M2.out)[1] <- colnames(a3.pairwise.out)[1]
> colnames(b3.M1.out)[1] <- colnames(a3.pairwise.out)[1]
> colnames(b3.M2.out)[1] <- colnames(a3.pairwise.out)[1]
> 
> 
> M <- nrow(a3.M1.out)
> J <- ncol(a1.M1.out)
> 
> rank.a3.M1 <- matrix(NA, M, J)
> rank.a3.M2 <- matrix(NA, M, J)
> rank.b3.M1 <- matrix(NA, M, J)
> rank.b3.M2 <- matrix(NA, M, J)
> rank.a3.pair <- matrix(NA, M, J)
> rank.b3.pair <- matrix(NA, M, J)
> 
> for (iter in 1:M){
+     rank.a3.M1[iter,] <- rank(a3.M1.out[iter,])
+     rank.a3.M2[iter,] <- rank(a3.M2.out[iter,])
+     rank.b3.M1[iter,] <- rank(b3.M1.out[iter,])
+     rank.b3.M2[iter,] <- rank(b3.M2.out[iter,])
+     rank.a3.pair[iter,] <- rank(a3.pairwise.out[iter,])
+     rank.b3.pair[iter,] <- rank(b3.pairwise.out[iter,])
+ }
> 
> 
> 
> rank.probs.a3.M1 <- matrix(NA, J, J)
> rank.probs.a3.M2 <- matrix(NA, J, J)
> rank.probs.b3.M1 <- matrix(NA, J, J)
> rank.probs.b3.M2 <- matrix(NA, J, J)
> rank.probs.a3.pair <- matrix(NA, J, J)
> rank.probs.b3.pair <- matrix(NA, J, J)
> 
> for (j1 in 1:J){
+     for (j2 in 1:J){
+         rank.probs.a3.M1[j1,j2] <- mean(rank.a3.M1[, j1] == j2)
+         rank.probs.a3.M2[j1,j2] <- mean(rank.a3.M2[, j1] == j2)
+         rank.probs.b3.M1[j1,j2] <- mean(rank.b3.M1[, j1] == j2)
+         rank.probs.b3.M2[j1,j2] <- mean(rank.b3.M2[, j1] == j2)
+         rank.probs.a3.pair[j1,j2] <- mean(rank.a3.pair[, j1] == j2)
+         rank.probs.b3.pair[j1,j2] <- mean(rank.b3.pair[, j1] == j2)
+     }
+ }
> 
> h.a3.M1 <- rep(NA, J)
> h.a3.M2 <- rep(NA, J)
> h.b3.M1 <- rep(NA, J)
> h.b3.M2 <- rep(NA, J)
> h.a3.pair <- rep(NA, J)
> h.b3.pair <- rep(NA, J)
> for (j in 1:J){
+     pvec <- rank.probs.a3.M1[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a3.M1[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.a3.M2[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a3.M2[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b3.M1[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b3.M1[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b3.M2[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b3.M2[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.a3.pair[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a3.pair[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b3.pair[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b3.pair[j] <- -1 * sum(holder)
+     
+ }
> 
> rank.a3.pair.hat <- apply(rank.a3.pair, 2, mean)
> ordvec <- order(rank.a3.pair.hat)
> 
> pdf(file="MM3-2-4-B-a3.pdf", height=6, width=7)
> plot(h.a3.M1[ordvec], xlab="Photo Index",
+      ylab="Entropy",
+      main="Black Coders", pch=0,
+      col="orangered1", ylim=c(0,5))
> abline(h=0)
> points(h.a3.M1[ordvec], pch=0, col="orangered1")
> points(h.a3.M2[ordvec], pch=1, col="dodgerblue")
> points(h.a3.pair[ordvec], pch=3, col="purple4")
> legend(x=0, y=5, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
+        col=c("orangered1", "dodgerblue", "purple4"))
> dev.off()
null device 
          1 
> 
> 
> 
> rank.b3.pair.hat <- apply(rank.b3.pair, 2, mean)
> ordvec <- order(rank.b3.pair.hat)
> 
> pdf(file="MM3-2-4-B-b3.pdf", height=6, width=7)
> plot(h.b3.M1[ordvec], xlab="Photo Index",
+      ylab="Entropy",
+      main="White Coders", pch=0,
+      col="orangered1", ylim=c(0,5))
> abline(h=0)
> points(h.b3.M1[ordvec], pch=0, col="orangered1")
> points(h.b3.M2[ordvec], pch=1, col="dodgerblue")
> points(h.b3.pair[ordvec], pch=3, col="purple4")
> legend(x=0, y=5, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
+        col=c("orangered1", "dodgerblue", "purple4"))
> dev.off()
null device 
          1 
> 
> 
> 
> ########################################################################
> ## END a3-b3 comparisons
> ########################################################################
> 
> 
> 
> 
> 
> ########################################################################
> ## START a4-b4 comparisons
> ########################################################################
> 
> ## load MCMC output
> load("MM3.2.a4.M1.out.Rda")
> a4.M1.out <- a4.M1.out[sample(1:nrow(a4.M1.out), nrow(a4.M1.out),
+                               replace=FALSE), ]
> load("MM3.2.a4.M2.out.Rda")
> a4.M2.out <- a4.M2.out[sample(1:nrow(a4.M2.out), nrow(a4.M2.out),
+                               replace=FALSE), ]
> load("MM3.2.a4.pairwise.out.Rda")
> a4.pairwise.out <- a4.pairwise.out[sample(1:nrow(a4.pairwise.out),
+                                           nrow(a4.pairwise.out),
+                               replace=FALSE), ]
> 
> load("MM3.2.b4.M1.out.Rda")
> b4.M1.out <- b4.M1.out[sample(1:nrow(b4.M1.out), nrow(b4.M1.out),
+                               replace=FALSE), ]
> 
> load("MM3.2.b4.M2.out.Rda")
> b4.M2.out <- b4.M2.out[sample(1:nrow(b4.M2.out), nrow(b4.M2.out),
+                               replace=FALSE), ]
> load("MM3.2.b4.pairwise.out.Rda")
> b4.pairwise.out <- b4.pairwise.out[sample(1:nrow(b4.pairwise.out),
+                                           nrow(b4.pairwise.out),
+                               replace=FALSE), ]
> 
> ## make MCMC output comparable
> a4.M1.out <- a4.M1.out[,1:48]
> a4.M2.out <- a4.M2.out[,1:48]
> b4.M1.out <- b4.M1.out[,1:48]
> b4.M2.out <- b4.M2.out[,1:48]
> 
> colnames(a4.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a4.M1.out))
> colnames(a4.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a4.M2.out))
> colnames(b4.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b4.M1.out))
> colnames(b4.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b4.M2.out))
> colnames(a4.pairwise.out) <- gsub("theta\\.", "", colnames(a4.pairwise.out))
> colnames(b4.pairwise.out) <- gsub("theta\\.", "", colnames(b4.pairwise.out))
> 
> for (j in 2:48){
+     a4.M1.out[,j] <- a4.M1.out[,1] + a4.M1.out[,j]
+     a4.M2.out[,j] <- a4.M2.out[,1] + a4.M2.out[,j]
+     b4.M1.out[,j] <- b4.M1.out[,1] + b4.M1.out[,j]
+     b4.M2.out[,j] <- b4.M2.out[,1] + b4.M2.out[,j]
+ }
> colnames(a4.M1.out)[1] <- colnames(a4.pairwise.out)[1]
> colnames(a4.M2.out)[1] <- colnames(a4.pairwise.out)[1]
> colnames(b4.M1.out)[1] <- colnames(a4.pairwise.out)[1]
> colnames(b4.M2.out)[1] <- colnames(a4.pairwise.out)[1]
> 
> 
> M <- nrow(a4.M1.out)
> J <- ncol(a1.M1.out)
> 
> rank.a4.M1 <- matrix(NA, M, J)
> rank.a4.M2 <- matrix(NA, M, J)
> rank.b4.M1 <- matrix(NA, M, J)
> rank.b4.M2 <- matrix(NA, M, J)
> rank.a4.pair <- matrix(NA, M, J)
> rank.b4.pair <- matrix(NA, M, J)
> 
> for (iter in 1:M){
+     rank.a4.M1[iter,] <- rank(a4.M1.out[iter,])
+     rank.a4.M2[iter,] <- rank(a4.M2.out[iter,])
+     rank.b4.M1[iter,] <- rank(b4.M1.out[iter,])
+     rank.b4.M2[iter,] <- rank(b4.M2.out[iter,])
+     rank.a4.pair[iter,] <- rank(a4.pairwise.out[iter,])
+     rank.b4.pair[iter,] <- rank(b4.pairwise.out[iter,])
+ }
> 
> 
> 
> rank.probs.a4.M1 <- matrix(NA, J, J)
> rank.probs.a4.M2 <- matrix(NA, J, J)
> rank.probs.b4.M1 <- matrix(NA, J, J)
> rank.probs.b4.M2 <- matrix(NA, J, J)
> rank.probs.a4.pair <- matrix(NA, J, J)
> rank.probs.b4.pair <- matrix(NA, J, J)
> 
> for (j1 in 1:J){
+     for (j2 in 1:J){
+         rank.probs.a4.M1[j1,j2] <- mean(rank.a4.M1[, j1] == j2)
+         rank.probs.a4.M2[j1,j2] <- mean(rank.a4.M2[, j1] == j2)
+         rank.probs.b4.M1[j1,j2] <- mean(rank.b4.M1[, j1] == j2)
+         rank.probs.b4.M2[j1,j2] <- mean(rank.b4.M2[, j1] == j2)
+         rank.probs.a4.pair[j1,j2] <- mean(rank.a4.pair[, j1] == j2)
+         rank.probs.b4.pair[j1,j2] <- mean(rank.b4.pair[, j1] == j2)
+     }
+ }
> 
> h.a4.M1 <- rep(NA, J)
> h.a4.M2 <- rep(NA, J)
> h.b4.M1 <- rep(NA, J)
> h.b4.M2 <- rep(NA, J)
> h.a4.pair <- rep(NA, J)
> h.b4.pair <- rep(NA, J)
> for (j in 1:J){
+     pvec <- rank.probs.a4.M1[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a4.M1[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.a4.M2[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a4.M2[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b4.M1[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b4.M1[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b4.M2[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b4.M2[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.a4.pair[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.a4.pair[j] <- -1 * sum(holder)
+ 
+     pvec <- rank.probs.b4.pair[j,]
+     holder <- pvec * log(pvec)
+     holder[pvec == 0] <- 0
+     h.b4.pair[j] <- -1 * sum(holder)
+     
+ }
> 
> rank.a4.pair.hat <- apply(rank.a4.pair, 2, mean)
> ordvec <- order(rank.a4.pair.hat)
> 
> pdf(file="MM3-2-4-B-a4.pdf", height=6, width=7)
> plot(h.a4.M1[ordvec], xlab="Photo Index",
+      ylab="Entropy",
+      main="Black Distractor Photos", pch=0,
+      col="orangered1", ylim=c(0,5))
> abline(h=0)
> points(h.a4.M1[ordvec], pch=0, col="orangered1")
> points(h.a4.M2[ordvec], pch=1, col="dodgerblue")
> points(h.a4.pair[ordvec], pch=3, col="purple4")
> legend(x=0, y=5, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
+        col=c("orangered1", "dodgerblue", "purple4"))
> dev.off()
null device 
          1 
> 
> 
> 
> rank.b4.pair.hat <- apply(rank.b4.pair, 2, mean)
> ordvec <- order(rank.b4.pair.hat)
> 
> pdf(file="MM3-2-4-B-b4.pdf", height=6, width=7)
> plot(h.b4.M1[ordvec], xlab="Photo Index",
+      ylab="Entropy",
+      main="White Distractor Photos", pch=0,
+      col="orangered1", ylim=c(0,5))
> abline(h=0)
> points(h.b4.M1[ordvec], pch=0, col="orangered1")
> points(h.b4.M2[ordvec], pch=1, col="dodgerblue")
> points(h.b4.pair[ordvec], pch=3, col="purple4")
> legend(x=0, y=5, legend=c("M1", "M2", "Pairwise"), pch=c(0,1,3),
+        col=c("orangered1", "dodgerblue", "purple4"))
> dev.off()
null device 
          1 
> 
> 
> 
> 
> ########################################################################
> ## END a4-b4 comparisons
> ########################################################################
> 
> 
> proc.time()
   user  system elapsed 
 17.559   1.949  19.514 
